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Abstract 

We  survey  the  advances  in  the  p-  and  the  h-p  versions  of  the  finite 
element  method.  An  up-to-date  list  of  references  related  to  these  methods  is 
provided. 
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1.  ^Introduction  and  brief  history. 

The  origins  of  the  finite  element  method  (FEM),  like  those  of  the 
spectral  method,  may  be  traced  back  a  long  time.  If  we  understand  the  FEM  as 
the  application  of  variational  principles  and  approximation  by  piecewise 
smooth  functions,  then  this  idea  was  already  used  by  Leibnitz  in  1696  (in 
one  dimension  with  piecewise  linear  functions).  In  two  dimensions.  Schellbach 


-Hrf  used  triangulation  and  piecewise  linear  functions, taee  also— [2]-Lq  Never- 

(V  f  h  ' 

theless,  the  modern  FEM  era  starts  with  the  paper  which  demonstrated  the 
potential  for  the  use  of  the  computer.  Since  then,  more  than  30,000  papers 


have  appeared (aee  [4tv  151,  [6] ) .(VThese  papers  are  generally  based  on  the 
^h-verslon  of  the  FEM,  where  the  accuracy  of  the  approximate  solution  is 
achieved  by  refining  the  mesh  while  using ^lo warder  polynomials  on  the  mesh 
The  spectral  method,  when  understood  to  be  the  use  of  variational 


.  ( 


principles  (or  other  methods,  such  as  collocation),  combined  with  the  use  of 
polynomials  of  high  degree,  was  already  known  to  Ritz.  A  method  of  this  type 
was  developed  (among  others)  by  Galerkln  [7]  and  was  discussed  in  detail  in 
[8],  for  eg..  A  number  of  further  developments  of  this  method,  attributed  to 
S.G.  Michlin,  may  be  found  in  his  many  papers  and  books.  For  example,  in  [9] 
he  discussed  principles  for  the  selection  of  basis  functions  and  outlined  a 
program  (based  on  polynomial  approximation)  for  the  Soviet  computer  M-20.  For 
various  theoretical  aspects  we  refer  to  the  important  paper  [10]. 

Both  methods  mentioned  were  found  to  have  their  strengths  and  weaknesses. 
The  FEM  provided  considerable  flexibility  and  was  well  suited  for  computer 
implementation.  The  spectral  method  offered  high  rates  of  convergence  when 
the  solution  was  smooth. 


In  the  1970s,  B. A.  Szabo,  recognizing  these  aspects,  suggested  and 
implemented  a  combination  of  the  two  approaches  to  utilize  the  advantages  of 
each.  Today,  this  combination  is  called  the  p-  and  the  h-p  version  of  the 
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finite  element  method.  If  the  mesh  Is  fixed  and  the  accuracy  of  the  solution 
is  achieved  only  by  increasing  the  degree  of  the  elements,  we  obtain  the 
p-verslon  of  the  FEM.  (If  the  domain  is  a  sruare  or  triangle,  and  is  under¬ 
stood  to  be  one  element,  then  the  p-version  is  identical  to  the  Ritz  method 
described,  for  exaunple,  in  [8]).  If  we  simultaneously  refine  the  mesh  and 
increase  the  degrees  of  elements  uniformly  or  selectively,  we  obtain  the 
h-p  version. 

The  first  theoretical  paper  addressing  the  p-version  was  [111  and  the  h-p 
version  [12].  Since  Szabo’s  original  work,  significant  progress  has  been  made 
for  these  methods  in  terms  of  theory,  implementation  and  engineering  applica¬ 
tions.  Some  of  these  achievements  are  addressed  in  this  paper. 

The  spectral- method  has  been  applied  extensively  in  the  last  15  years  to 
problems  in  fluid  mechanics.  Recently,  there  has  been  interest  shown  in  using 
this  method  over  partitioned  domains  (rather  than  a  single  one,  see  for  eg. 
[131).  In  this  context,  the  spectral  method  over  a  partitioned  domain  is  very 
similar  to  the  h-p  version,  though  the  emphasis  of  the  two  methods  is 
different  -  the  h-p  version  of  the  FEM  concentrating  on  the  special  needs  of 
structural  mechanics  analysis,  while  the  spectral  method  being  specialized 
more  for  fluid  mechanics. 

There  are  many  prograuns  based  on  the  h-version  of  the  FEM,  some  major 
commercial  ones  being  MCNASTRAN,  ADINA,  ANSYS  and  others.  There  are  only  two 
commercial  programs  based  on  the  p-  and  h-p  version,  FIESTA  and  PROBE,  in 
addition  to  a  large  research  program  called  STRIPE.  Other  commercial  prograuns 
based  on  the  p-version  are  being  developed  at  various  places  and  will  be  on 
the  market  in  the  near  future.  The  authors  of  this  paper  have  experience  with 
PROBE  and  references  to  it  (rather  thaui  ary  alternative)  are  for  convenience 


only. 
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The  h-p  version  of  the  finite  element  method  has  various  features  which 
are  reflected  in  the  implementation  and  architecture  of  the  program  and  are 
different  from  the  h-version.  Recent  advances  in  computer  hardware  bring  to 
the  forefront  the  problem  of  reliability  of  computations  and  the  repercussions 
of  the  rapidly  changing  ratio  between  computer  and  human  costs.  The  h-p 
version  offers  various  essentially  new  possibilities  which  the  h-version  does 
not.  As  examples,  we  mention  the  (a-posteriori )  assessment  of  the  errors  of 
the  FEM  calculations,  new  possibilities  in  the  modeling  of  plates  and  shells 
(with  models  of  Reissner-Mindlin  type  being  naturally  created)  and  inherent 


parallelization.  In  addition,  the  h-p  version  shows  remarkable  robustness, 
for  eg.  with  respect  to  locking  phenomena. 

"-''-'4  In  this  paper  we^present^i  survey  of  the  state  of  the  art  of  the  p  and 
h-p  versions.  The  emphasis  is  on  the  theoretical  aspects  related  to  their  use 
in  approximating  elliptic  equations  stemming  from  structural  mechanics. 


2.  The  model  problem. 

Problems  in  structural  mechanics  and  the  mechanics  of  solids  are 
typically  characterized  by  elliptic  partial  differential  equations  with 
piecewise  analytic  data,  pertaining  to  the  boundary,  boundary  conditions, 
coefficients  and  right  hand  side.  Consequently,  one  can  expect  special 
features  in  the  solution  which  should  be  somehow  exploited  by  the  numerical 
method  used.  In  this  section,  we  mention  some  typical  available  results.  For 
simplicity  and  brevity,  we  restrict  our  discussion  to  the  two  dimensional 
case. 

2 

Let  Qj  c F  ,  J  ■  1,2,  •••,M  be  simply  connected  domains  with  boundaries 

B  j 

8Q.  =  r(J)  *  U  fj^.  f  are  analytic  simple  arcs  which  we  call  edges, 

J  1-1  1 

while  the  ends  of  the  edges  are  called  vertices.  If  nij  *  3, 


3 


« 


respectively  =  4,  we  call  flj  a  curvilinear  triangle,  respectively 

quadrilateral ■  Otherwise  It  will  be  a  curvilinear  polygon.  We  will  assume 

that  O.nO.  is  either  empty,  is  a  common  vertex,  or  is  a  common  edge.  Let 
J  M 

n  be  the  interior  of  U  Q.  and  assume  that  Q  is  a  domain.  By  r  we 

J-l  J 


denote  the  boundary  of  Q. 


The  edges 


not  belonging  to 


r°  will  be 


called  Interface  edges . 

We  now  consider  a  model  problem  for  second  order  scalar  elliptic 
differential  equations  written  in  the  weak  form.  Let 

(2.1)  B(u,v>  -  J  (  £  j  Ssr  1^7  ] 

n  i , j=i, 2  J 

where  a^  =  a^  are  analytic  functions  on  0^  satisfying  the  standard 
elliptlclty  condition 


(2  2)  [ 

l.J-1.2 

Further,  let 

(2.3)  Ft(v)  *  J  fv  dx^ 

n 

where  f  Is  an  analytic  function  on  n^,  l  ■  1,***,M.  Let  u  be  continuous 
on  D  ■  U  F. ,  F  »  F^cr0  and  analytic  on  F..  D  will  be 

r  jSD  1  J 

called  the  Dlrichlet  boundary. 

Finally,  let  H  *  r^-D  be  the  Neumann  boundary  and  let  g  be  defined 
on  H  and  analytic  on  every  1^  c  H,  with 

(2.4)  f2(v)  ■  £  J  g  v  ds. 

ri‘K  ri 

The  exact  solution  of  our  model  problem  is  defined  in  the  usual  way:  find 
uQ e H*(Q),  Uq  ■  u  on  D  such  that 
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(2.5) 


B(uQ,v)  =  (F1  +  F2)(v) 

for  all  ve°H*  =  {veH*(fl),  v  =  0  on  D> .  If  D  Is  empty,  then  the  usual 

l/? 

solvability  condition  for  has  to  be  satisfied.  By  fluHg  ■  (B(u,u)) 

we  denote  the  energy  norm  of  u.  It  is  equivalent  to  the  H*(Q)  seminorm 
(and  to  the  H1^!)  norm  on  °H1  (Q) ) .  The  most  important  cases  are  when  a^ 
are  constant  on  every  l  -  1 ,  •  •  • ,  M. 

A  similar  formulation  holds  for  the  elasticity  problem.  Here,  u  =  ( u ^ , u_ ) 

and 

C2.6)  B(u.v)  =  J  [  [  bijklc1j(-*€kl^-^dXldX2 

n  i,J,k, 1=1,2 

cu(s>  -  i 

with  the  standard  assumptions  about  bjj  and  the  functional  F.  The 
assumptions  about  the  piecewise  analyticity  are  analogously  formulated  as 
before.  Here  the  boundary  conditions  can  be  more  general,  combining  traction 
components  and  displacements  on  particular  Tj. 

In  the  case  of  an  isotropic  material,  the  bilinear  form  is 

(2.7)  B(u.v)  *  J  [  l  e1J(u)eij(v)  ♦  ^  (dlv  u  div  v)  jdx^. 

0  i.J-1,2 

Here,  v  is  the  Poisson  ratio,  0 < v  < |  and  E  is  the  Young's  modulus 
of  elasticity.  The  form  degenerates  for  v  — but  regularity  of  the 
solution  is  preserved.  A  similar  degeneration  with  preservation  of  the 
regularity  properties  also  occurs  for  general  anisotropic  materials  (see  [14], 
[ 15] ,  for  eg. ) . 

The  behavior  of  the  solution  is  essentially  very  similar  in  both  the 
scalar  case  and  the  elasticity  problem.  The  solution  u^ 
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a)  is  analytic  in  Q,-  U  aJ*^,  i  *  1, •••,*,,  J  *  1,***,M. 

J  i.j  J 

b)  has  a  special  singular  behavior  in  the  neighborhood  of  every  A^. 

The  behavior  of  the  solution  is  best  understood  when  the  operator  in  our 
problem  is  the  Laplacian.  In  this  case,  near  a  vertex,  we  have 

J  S  N 

a  .+m 

(2'8)  U0  ■  I  l  l  CJS»  *js.(8,r  108  r  *  V 

J=1  s=0  m=0 

Here,  (r,6)  are  local  polar  coordinates  at  the  corner  point  under  considera¬ 
tion.  The  decomposition  (2.8)  is  such  that  the  remainder  u^  is  smoother 
than  the  terms  in  the  sum.  The  functions  ^jgm  sure  smooth  (in  our  case 
piecewise  analytic).  We  have  S  *  0  except  for  special  cases,  when  S  *  1 
is  possible.  N  may  be  0  or  positive.  For  eg.,  N*0  in  the  cases  when 
a^  are  nonconstant  or  T*^  are  curved.  For  details  we  refer  to  [16J.  The 
norm  of  uQ  in  (2.8)  depends  on  the  geometry  and  diverges  to  »  when  the 
geometry  converges  to  certain  exceptional  cases.  The  coefficients  cjsm  are 
related  to  the  stress  intensity  factors.  They  can  be  global,  depending  on 
Uq,  or  local,  depending  only  on  the  input  data  at  the  vertex. 

In  the  case  of  the  elasticity  equations,  the  results  are  similar, 
although  not  as  detailed.  The  coefficients  in  (2.8)  can  now  be  complex 

and  the  conditions  for  S  and  M  are  not  completely  characterized. 

There  is  often  a  practical  need  in  actual  problems  to  know  the  values  of 

a.  and  the  functions  i/i .  In  [17]  a  general,  adaptive,  completely  robust 
J  jsm 

algorithm  for  determining  these  coefficients  and  functions  is  given.  As  an 
example  of  the  complicated  structure  of  these  coefficients,  we  show  the  two- 
material  (anisotropic)  case  when  zero  displacements  are  prescribed  at  the 
boundary  (Fig. 2.1). 
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ADHESIVE 


GRAPHITE 


Fig. 2.1  the  scheme  of  the  two  material  domain. 


The  two  materials  are  typical  anisotropic  ones  used  In  engineering. 

Graphite  is  highly  anisotropic  while  adhesive  Is  only  slightly  anisotropic. 

Figs.  2.2  show  the  first  six  (or  seven)  with  smallest  real  part  as 

-5 

functions  of  the  angle  (the  accuracy  is  10  ).  Fig.  2a  shows  the  real  part 

of  a..  If  a  Is  complex,  values  are  denoted  by  circles.  Fig.  2b  depicts 
J  J 

the  imaginary  part.  For  details  see  [17]. 


0  20  40  40  N  101  110  14  IfO  110 

ANGLE  OMEGA]  Ot  DCOMB 


Fig. 2. 2a)  The  real  part  of  a, 


»  40  60 


120  MO  160  1)0 


ANCLE  OMEGA!  IN  I 


Fig. 2. 2b)  The  Imaginary  part  of  a.. 


There  Is  a  vast  literature  devoted  to  the  analysis  of  the  decomposition 
(2.8).  We  mention  here  [18],  [19],  [20],  [21]  and  references  given  therein. 

The  above  decomposition  is  valid  when  the  input  data  are  smooth  but  not 
necessarily  analytic.  Another  characterization  of  the  regularity,  given  In 
terms  of  countably  normed  spaces,  was  analyzed  In  [22],  [23],  [24],  [25].  For 
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example,  It  Is  shown  that  for  the  exact  solution  u^. 


(2.9) 


I  D“u0 1 


SC 


£pl«l+0-lj 


-1 


'a! 


0  <  0  <  1 


holds  for  all  ail.  Here  C  and  d  are  Independent  of  a.  In  the  above 
references,  the  complete  characterization  of  these  normed  spaces  together  with 
corresponding  trace  spaces  and  extensions  is  analyzed. 

In  3  dimensions,  the  situation  is  more  complicated  due  to  the  presence  of 
both  edge  and  vertex  singularities.  As  a  result,  the  regularity  may  be 
different  along  different  directions,  leading  to  the  use  of  anisotropic 
spaces. 

Let  us  finally  point  out  that  the  regularity  of  the  solution  of  our  model 
problem  may  also  be  characterized  in  terms  of  standard  Sobolev  (or  Besov) 
spaces.  Accordingly,  if  a  is  the  minimum  over  all  vertices  Aj  of  the 
exponents  (or  ReCa^)  if  aj  are  complex)  in  (2.8),  then  we  have 

(2.  10)  uchNq)  V  k  <  1  +  a. 

Most  classical  finite  element  error  estimates  rely  on  regularity  results 
of  the  form  (2.10). 


3.  The  h-,  p-  and  h-p  version  of  the  finite  element  method. 

To  illustrate  the  basic  results,  we  restrict  ourself  to  very  special 
canes,  although  the  available  results  are  completely  general. 

Let  vis  consider  the  case  when  n  is  an  L-shaped  domain,  as  shown  in 
Fig. 3. 1. 
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B  A 


Fig. 3.1  The  scheme  of  the  L-shaped  domain. 

We  consider  the  elasticity  problem  (2.7)  with  f  =  0  and  traction 
boundary  conditions  such  that 

(3.1)  U1  “  iG  r*(  (k~Q(X+1))  cos  X0-X  cos  (X— 2)©] 

U2  =  IG  (k+q(*+1^  sin  +  *  sin  ( X— 2 )© ] 

where 

k  *  3  -  4i>, 

G  *  E/2( 1+v) 

„  _  _  X  sin  (X-l)«/2 
W  "  sin  (X+Uu/2 

with  v  *  0.3,  X  =  0 . 544484 ,  «  =  |  k.  The  sides  OE,  and  OA  (see  Fig. 3.1)  are 
traction  free.  Solution  (3.1)  is  one  term  in  the  decomposition  mentioned  in 
Section  2. 

As  usual,  we  introduce  a  mesh  to  partition  Q.  For  simplicity,  we  first 
consider  the  case  of  a  uniform  partition,  characterized  by  the  parameter  h 
(see  Fig. 3. 2). 
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The  finite  element  spaces  VcH^(U)  will  consist  of  continuous  piecewise 
polynomials  of  degree  p  on  the  squares  of  the  mesh.  The  exact  set  of 
polynomials  used  over  each  square  consists  of  either  Q^,  the  set  of 
polynomials  of  degree  p  separately  in  each  variable,  or  Q^,  the  minimal  set 
containing  polynomials  of  total  degree  p.  See  eg. [26]  for  details.  The 
finite  element  solution  u^  is  then  defined  as  usual  by 

Btu^.v)  -  F(v)  Vve  V 

with  the  error  satisfying 

i,h0-Hfeiie  *  BiJ  "V^E- 

WcV 

The  space  V  (and  hence  also  Up^)  Is  characterized  by  two  parameters,  p 
and  h,  so  that  V  *  V(p,h).  By  N(p, h)  we  denote  the  dimension  (l.e.  the 
number  of  degrees  of  freedom)  of  V(p,h).  In  order  to  obtain  a  desired 
accuracy  for  our  approximation,  we  use  am  extension  procedure,  l.e.  a 
procedure  to  Increase  the  dimension  N(p,h).  This  can  be  of  three  types: 

a)  h-verslon.  Here  p  is  fixed,  usually  at  a  low  level  (eg.  1  or  2)  and 
we  achieve  the  desired  accuracy  by  taking  h— >0. 

b)  p- version.  In  this  case  h  Is  fixed,  l.e.  the  same  mesh  Is  used  and 
p— »«,  l.e.  the  accuracy  Is  achieved  by  increasing  p. 
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c)  h-D  version.  Here  h  and  p  are  simultaneously  changed  (either 
uniformly  or  selectively). 

In  general,  we  will  be  Interested  In  the  relative  error  ||e||_  « 

ER 

||uQ  -  Upgll^lluQllg.  We  present  some  theorems  for  our  approximation  which  are 
the  special  cases  of  those  proven  In  [27],  [28]. 

Theorem  3. 1  [27] 

(3.2)  ||u0-uFEllESChminU*p)p"2X 

where  C  Is  a  constant  Independent  of  h  and  p. 

The  above  theorem  holds  for  both  choices  Q  and  Q' .  In  2  dimensions, 

P  P 

-2  2 

N  *  h  p  so  that  to  obtain  the  optimal  asymptotic  rate  minimizing  N(p, h),  we 
choose  h  =  1,  i.e.  the  p-verslon.  (Of  course,  N  does  not  completely  measure 
the  needed  work.  Moreover,  the  accuracy  measured  In  the  energy  norm  Is  not 
necessarily  the  accuracy  we  are  seeking  In  practice).  Figs.  3.3a  and  3.3b 
show  the  errors  for  the  h-  and  p-  versions  respectively  with  elements  of 
type. 
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The  figures  are  drawn  in  the  log  log  scale.  They  show  that  (3.2)  correctly 
characterizes  the  asymptotic  behavior  which  is  defined  by  the  slope  in  the 
figure. 

Let  us  now  consider  a  non-uniform  geometric  mesh  with  n  layers  and  with 
ratio  0.15.  This  is  shown  in  Fig. 3. 4  for  n  *  2. 


detail 
n*2 

Fig. 3. 4  The  geometric  mesh 

The  shape  functions  are  now  the  usual  mapped  polynomials  using  blending 
mapping  techniques  for  the  circular  sides.  For  details,  see  eg.  [29]. 

In  Fig. 3. 5  we  show  the  error  for  different  numbers  of  layers  and  the 
(uniform)  degree  p  in  the  log  log  scale 


,2 


I* 


Fig. 3. 5  The  error  lleSER  as  a  function  of  n  and  p. 

1/3 

Fig. 3. 6  shows  the  error  behavior  In  log  ||e||^xN  for  selected 
combinations  of  n  and  p. 


Fig. 3. 6  The  error  ll®#^  Tor  selected  combinations  of  (n,p). 

-T3vfi 

We  see  that  HeH^  *  Ce  Estimates  of  this  type  have  been  analyzed 

in  [30],  [31],  [32],  employing  the  smoothness  characterization  of  (2.9).  We 
have 

Theorem  3.2  [32].  Let  the  mesh  with  n  layers  be  considered  and  let 

pn  S  p  £  un,  p  >  1  0  <  p  <  v  <  «. 

Then  If  uQ  satisfies  (2.9),  we  have 

0yo-yFEllESCe"r  ^  •  r>0 

In  actual  computations,  the  remeshing  required  by  the  h-p  version  is 
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a  disadvantage  due  to  increased  human  cost.  Consequently,  the  usual  practice 
while  employing  p  and  h-p  versions  codes  is  to  use  a  fixed,  strongly 
refined  mesh  and  then  increase  p  (i.e.  the  p-version).  The  mesh  design 
should  be  such  that  the  desired  accuracy  will,  as  far  as  possible,  be  achieved 
for  the  optimal  pair  (n, p).  In  [33],  [34]  an  attempt  to  design  an  expert 
system  for  such  selection  is  presented. 

Theorem  3.1  is  a  special  case  of  the  following,  more  general  theorem, 
proved  in  [27]. 

Theorem  3.3  [27]  Assume  that  u^  has  the  form  (2.8)  with  a  =  min  a^.  Then 

r  min(a,p-ak 

H Uq-UfeII e  *  Cg(h'P>s)  min  h  > - 2a -  ’  8(h,p,S)  -  max  ( |  log  hi  ,  I  log  pi  ) 

l  p  J 

Theorems  3.1  and  3.3  show  the  interesting  fact  that  the  convergence  rate 
of  the  p-verslon  is  twice  that  of  the  h-version  when  a  uniform  or  (more 
generally)  a  quasiuniform  mesh  is  used.  This  fact  was  proved  in  [11]  for  the 
p-version  (see  also  [35]). 

In  this  connection,  the  following  result  from  [27]  is  useful  when  u^  is 
only  known  to  be  in  some  Sobolev  space  as  in  (2.10). 

Theorem  3.4  [27].  Let  UqCH^UI),  k>  1.  Then  if  the  spaces  V  * 

V(p, h)  are  based  on  a  uniform  (or  quasiuniform)  family  of  meshes, 

0.3)  IS0-SFElE»ChV(lt'I)l!i0lHk(n) 

where  p  -  min  (p,  k-1)  and  C  is  independent  of  uQ,h,p. 

Theorem  3.4  Improves  the  classical  estimate 

l,^0"iWlE*C(p)hM 

for  the  h-verslon  by  explicitly  showing  how  the  constant  C(p)  decreases  when 
p  is  increased.  Note  that  theorem  3.3  Is  a  more  refined  result  for  solutions 
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of  the  form  (2.8),  since  using  (2.10)  with  theorem  3.4  will  not  yield  the 
observed  doubling  in  the  rate  of  convergence  of  the  p-version. 

The  p-version  has  been  analyzed  for  3-dimensional  problems  in  (36),  [37]. 

Various  problems  arising,  for  example,  from  the  theory  of  plates  and 
shells  may  be  described  by  elliptic  equations  of  order  2a  where  m> 1.  For 
such  problems,  if  the  elements  used  are  conforming  piecewise  polynomials,  then 
they  must  have  m-1  continuous  derivatives  over  Q.  Approximation  results 
for  the  p-verslon  using  such  Cm  1  elements  have  been  established  in  [38], 
where  it  is  shown  that  once  again,  due  to  the  presence  of  ra  type  singulari¬ 
ties  in  the  solution,  one  obtains  twice  the  rate  of  convergence  of  the 
h-version.  The  case  m=2  was  originally  discussed  in  [39]  where  some 
computational  results  using  C1  elements  are  presented.  Theorems  for  the  h-p 
version  for  equations  of  order  2a  are  given  in  [40]. 

For  second  order  problems,  the  results  in  theorems  3. 1-3.4  hold  not  only 
for  square  or  triangular  elements  but  also  for  curvilinear  elements  having 
some  uniformity  properties  with  respect  to  their  mapping  onto  standard 
elements.  For  details,  see  eg.  [31]. 

As  we  have  seen,  the  presence  of  singularities  significantly  decreases 
the  rate  of  convergence.  In  addition,  the  p-version  is  influenced  by  the 
"pollution"  problem  [41].  By  this  we  mean  the  effect  of  an  error  in  one 
element  (usually  due  to  a  singularity  present  in  the  true  solution  over  that 
element)  permeating  into  adjacent  elements  (where  the  exact  solution  is 
regular).  The  pollution  problem  is  more  serious  if  the  stresses  are  of 
interest.  It  may  be  overcome  by  using  refined  meshes  (a  few  layers)  in  the 
area  of  singularity.  Another  approach  to  deal  with  this  problem  is  to  use 
properly  mapped  shape  functions.  For  details,  see  [42]. 

So  far,  we  have  assumed  in  our  model  problem  that  the  Dirlchlet  boundary 
set  is  empty.  The  theorems  we  mentioned  above  are  valid  without  changes  when 
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u  *  0  on  D  i.e.  the  Dirichlet  conditions  are  homogeneous.  Then  we  simply 
use  °V(p,h)  «  V(p,h)  instead  of  V(p,h). 

In  the  case  of  nonhomogenous  essential  boundary  conditions,  we  have  to 
approximate  u  by  Wpg  so  that  «pg  is  in  the  trace  space  of  the  finite 
element  space  V(p, h).  This  is  done  by  a  projection  in  the  Hy ( r ^ )  norm, 

0  57  SI.  The  cases  7  *  1,1/2  have  been  analyzed  in  [28],  [43]  and  the 
general  case  in  [44].  The  results  show  that  for  smooth  w,  the  optimal  rate 
of  convergence  is  achieved  when  l/2SySl.  More  precisely.  Theorem  3.4 
holds  for  I/2S7SI  and  wcH^r-j),  s>l.  For  h>  unsmooth,  e.g. 
ueH^Tj),  |  <  s  <  |,  the  optimal  rate  of  convergence  using  the  projection 
approach  has  been  established  only  for  7  ■  |.  Numerical  results  cure  given  in 
[30]  and  [44].  This  problem  does  not  occur  with  the  h-p  version  when  u  is 
singular  in  the  neighborhood  of  the  vertices.  In  [45],  we  have  analyzed  a 
class  of  constrained  boundary  conditions  which  are  important  in  practice  in 
structural  mechanics. 

So  far  we  have  only  mentioned  the  solution  of  elliptic  problems.  The  p 
and  h-p  version  can  also  be  used  for  other  types  of  problems.  For  example. 
In  [46],  [47]  we  analyze  the  method  for  solving  parabolic  equations  when  the 
h-p  version  is  used  in  both  the  time  and  space  variables. 

4.  The  problem  of  optimal  meshes  and  adaptive  approaches. 

In  the  previous  section,  Fig.  3.3a  showed  the  convergence  rate  for  the 
h-verslon  using  a  uniform  mesh.  This  rate  may  be  improved  by  using  better 
meshes  in  certain  cases. 

The  problem  of  optimal  meshes  for  the  h-  and  h-p  versions  was  studied  in 
detail  for  1  dimension  In  [48]  and  for  2  dimensions  in  [49].  Let  us  mention 
some  one  dimensional  results. 

We  consider  the  simple  model  problem 
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-u"  «  f,  u(0)  ■  u(l)  ■  0 


ot  1 

with  the  exact  solution  u(x)  -  x  -  x,  a>5-  Let  x,  denote  the  mesh 
points.  For  the  h-version,  the  radical  mesh  x^  ■  -  i  ■  l,--*,m  is 

optimal. 

p  ♦  i 

Theorem  4.1  [48].  The  radical  mesh  with  0  *  - -  is  optimal  and 


1 


lim  mpl|u-u_||E  *  C( 
m-w  *■ 


fp  *  5] 

*.  p)-* — r1 


p+5 


1 

a  i 


where 


C(a,p) 


aT(a) Isin  wa|  r(p-a+l) 
V*  4P/2p+l  T( p+1/2) 


Theorem  4. 1  shows  that  for  the  h-version,  the  best  possible  rate  of 
convergence  is  0(hp),  which  is  algebraic,  and  not  exponential  (as  for  the  h-p 
version) . 

For  the  h-p  version,  the  optimal  mesh  Is  a  geometric  one  with  ratio  q, 

uj—  i 

xt=q  ,  0<q<l,  1  ■  1,2,  •••,m  and  the  optimal  p-dlstrlbutlon  is 
linear,  p1  ■  [si]  +  1  where  [a]  denotes  the  Integral  part  of  a  and 
Pt  is  the  degree  of  the  element  (x^.x^).  s  will  be  called  the  slope. 

In  the  case  of  uniform  p,  we  use  p  «  [sm]  +  1.  The  optimality  is 
understood  in  the  sense  that  the  error  using  the  optimal  mesh  and  optimal 


Vlft*l/2  Jii 

degree  distribution  has  the  same  exponential  rate  qQpt  as  the  best 


achievable  rate  among  all  mesh  and  p  distributions  with  the 
degrees  of  freedom  N. 


number  of 


Theorem  4.2.  We  have  qQpt  -  (VS-lr  *  0.17  and  sQpt  ■  2a- 1.  Then 


'“o-'Ve'  s 
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In  the  case  of  uniform  degree  distribution  we  have 


Theorem  4.3. 


%pt 


(VS-l)2,  sQpt  *  2a- 1  and 

'“o->ve'e  s  c'*><‘r,/4'N  vr°"2 


<r  ■  min  (2a-  l,a). 


We  have  seen  that  for  p  uniform,  the  radical  mesh  is  optimal.  Hence  we  can 

ask  about  the  envelope  of  optimal  radical  meshes.  Then  the  radical  meshes 

-4/e2 

tend  to  a  geometric  one  with  ratio  q  »  e  «*0.54  and 

2 

s  =  4(a-l/2)/e  *  0.54(a-l/2).  These,  together  with  many  more  detailed 

results  as  well  as  numerical  experimentation  are  given  in  [48]. 

In  two  dimensions,  the  situation  is  more  complicated.  Nevertheless,  the 
linear  distribution  of  p  and  geometric  mesh  are  once  again  optimal.  The 


— ~tVN 

estimates  will  be  of  the  type  e  in  contrast  to  e  *  in  one  dimension. 

For  circular  elements  and  optimal  choice  of  the  degrees  of  elements  in 
different  directions,  we  can  achieve  the  rate  e  *  too.  For  details  we 
refer  to  [49], 

The  above  results  indicate  that  the  geometric  mesh  with  q  *  0. 15  is  the 
right  mesh  for  practical  use.  It  is  preferable  to  over-refine  the  mesh 
slightly.  The  selection  of  the  number  of  layers  can  be  made  in  an  expert 
system  mode  or  adaptively.  Adaptive  approaches  were  addressed  in  [49],  [50]. 
Let  us  mention  that  the  codes  FIESTA  and  STRIPE  have  some  adaptive  features 
with  shape  functions  being  selected  in  an  anisotropic  way.  In  PROBE,  the 
determination  of  the  p-distribution  is  done  at  present  by  the  user. 


5.  The  p  and  h-p  version  for  integral  equations  and  mixed  methods. 

The  theorems  in  Section  3  (and  4}  were  based  on  approximation  theory 
results  in  the  H1  norm.  We  can  proceed  analogously  for  cases  where  we  have  a 
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coercive  bilinear  form  over  some  other  space  H  for  which  corresponding 
regularity  and  p-  and  h-p  version  approximation  results  are  known.  In  [511 
this  procedure  is  extended  to  the  boundary  element  method.  Consider  for 
example  the  model  problem  from  Section  2,  of  Laplace’s  equation  on  a  polygonal 
domain  when  both  the  Oirichlet  and  Neumann  sets  are  present.  Then  the  problem 
cam  be  formulated  on  T*1  (see  [51])  as  a  system  of  integral  equations  with  the 
unknowns  being  given  by  the  pairs  ,  |^lj-  j  where  F^  c  tt,  c  D.  This 

problem  may  be  put  in  the  form  B(u, v)  ■  F(v).  It  satisfies  a  Carding 
inequality  which  is  sufficient  to  obtain  am  optimal  rate  of  convergence  for 
the  Galerkln  approximation.  In  [52],  the  h  and  h-p  versions  were  analyzed  and 
it  warn  shown  that  the  rate  of  convergence  of  the  p-verslon  is  twice  as  high  as 
that  for  the  h-ve'rslon  (with  uniform  mesh),  similar  to  the  cases  discussed 
earlier.  We  cam  also  obtain  an  exponential  rate  of -convergence  for  the  h-p 
version  with  a  properly  chosen  (geometric)  mesh  and  degrees  amalogously 
selected  as  In  the  previous  sections.  For  details,  see  [53],  [54],  [55].  For 
adaptive  procedures  in  the  h-p  version  for  Integral  equations,  we  refer  to 
[56]. 

The  problems  discussed  so  far  have  been  stable  (in  the  sense  that  they 
are  coercive  or  a  Carding  inequality  holds).  For  mixed  methods,  one  must 
first  establish  the  stability  of  the  approximate  subspaces  used,  via  an 
inf-sup  condition.  The  stability  of  the  p-version  in  the  context  of  certain 
mixed  methods  for  Stokes’  problem  has  been  discussed  in  [57],  [58].  (See  also 
spectral  method  references).  In  [59],  the  Raviart-Thomas  and  the  Brezzi- 
Douglas-Marlni  spaces  for  the  mixed  formulation  of  linear  elliptic  problems 
have  been  shown  to  be  stable  and  possess  optimal  convergence  properties  in 
terms  of  the  h-p  extension  using  quasi uni form  meshes.  The  p-extenslon  of 
Raviart-Thomas  elements  for  quasi linear  problems  is  analyzed  in  [60]. 
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6.  The  h-p  version  and  mathematical  Modeling. 


Let  vis  consider  as  a  model  problem  the  problem  of  plate  bending.  This 

problem  Is,  in  fact,  a  3  dimensional  problem  over  a  "thin"  domain 

3 

12  *  wx  (-t/2,  t/2)  c  R  .  Two  dimensional  formulations  such  as  Kirchhoff  or 
Relssner-Mindlin  models  (among  others)  are  dimensionally  reduced  formulations 
of  this  problem.  These  formulations  generally  are  asymptotically  Identical 
for  t  — »0,  but  yield  different  results  for  t > 0.  The  solutions  of  these  2 
dimensional  models  have  to  be  understood  as  approximations  of  the  3  dimen¬ 
sional  formulation.  The  error  depends  on  the  type  of  input  data  (for  eg., 
clamped  or  simply  supported  plate),  thickness  and  the  aim  of  the  computation. 
The  h-p  version  gives  a  natural  tool  which  leads  to  a  hierarchical  set  of 
formulations.  Denoting  the  displacements  by  u  *  (u^, 
reduction  can  be  understood  as  a  projection  on  the  space  of  solutions  of  the 
form 


Ug.u^),  the  dimensional 


(6.  1) 


s. 


U1*X1,X2,X3*  “  E  ul^  (xl*x2)x3  ^ 


J-l 

s„ 


U2(X1,X2,X3^  =  E  ^  ^xl,x2^x3  ^ 


J-l 


WVV  "  E  U3J)(X1’X2)X32J  ' 

J-o 


For  example  If  v  ■  0  (y  »  Poisson  ratio)  then  the  choice  Sj  ■  Sg  ■  1,  s3  **  0 
leads  to  the  Relssner-Mindlin  model.  For  y>0  one  has  to  take  *  1.  For 
more  about  dimensional  reduction  we  refer  to  [61]  and  references  therein.  The 
error  of  the  reduced  formulations  depends  on  various  factors.  For  example 
(see  [62]),  for  the  simply  (soft)  supported  uniformly  loaded  Relssner-Mindlin 
plate  with  angle  30*.  side  length  1  and  thickness  t  *  0*1,  0*01,  the  errors 
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in  the  energy  norm  are  4.34%  and  15.41%  respectively.  The  form  (6.1)  can  be 
understood  as  the  p-version  with  the  polynomial  degrees  in  the  x^  direction 
being  different  from  those  in  x^.Xg.  From  this  point  of  view,  the  h-p 
version  is  a  natural  tool  for  deriving  plate  models  and  assessing  their  error 
(see  next  section).  The  program  PROBE  has  these  types  of  features  for  appli¬ 
cation  to  plates  and  shells,  as  well  as  for  transitions  where  s  is  changed 
in  various  parts  of  the  domain.  As  shown  in  [62]  and  [63],  the  various 
boundary  conditions  (hard,  soft)  have  a  significant  Influence  on  the  solution. 
For  more,  we  also  refer  to  [64]. 

7.  Extraction  techniques. 

Usually  in  computational  practice,  the  solution  u  of  our  variational 
problem  is  only  a  tool  to  get  the  primary  quantity  of  interest.  For  example, 
the  goal  of  the  computation  may  be  to  find  the  stresses  at  a  point,  or  the 
maximal  stress  (e.g.  Mlses  equivalent  stress)  over  a  region,  or  the 
resultants  (reactions,  moments)  in  the  plates  and  shells,  stress  intensity 
factors,  etc.  Mathematically,  we  are  Interested  in  evaluating  the  values  of 
certain  functionals.  This  can  be  done  in  a  trivial  way  (for  example,  by 
differentiating  the  finite  element  solution  u^g)  or  using  more  sophisticated 
approaches  which  lead  to  more  accurate  results  (with  accuracy  being  of  the 
order  of  the  error  in  the  energy,  rather  than  the  energy  norm).  Such 
techniques,  called  extraction  techniques,  were  addressed,  for  e.g.,  in  [65] 
and  applied  in  various  Important  contexts  (see  [66],  [67]). 

As  an  example  of  an  extraction  technique,  we  present  computations  for  the 
stress  intensity  factors  for  the  model  problem  Introduced  in  Section  3,  but 
with  the  exact  solution  consisting  of  two  terms  of  the  expansion,  Aj  ■  0.54448 
and  *2  ■  0.90853.  We  have  selected  the  intensity  factors  otj  ■  1  (mode  1) 


and  ctg  ■  2  (mode  2).  Fig. 7.1  shows  the  error  of  and  a ^  in  the  scale 

1/3 

log  exN  for  the  two  layer  mesh  as  well  as  the  error  In  the  strain  energy 
(not  energy  norm). 


POLYNOMIAL  DEGREE 


Fig. 7.1  Convergence  of  the  stress  intensity  factors  computed  by  the 
extraction  technique. 

We  see  that.  In  fact,  the  accuracy  in  the  stress  Intensity  factors  Is  of 
the  same  order  as  that  In  the  strain  energy  and  that  for  high  p,  the  second 
mode  Is  more  accurate  than  the  strain  energy  (see  [65]  for  an  analysis).  As 
we  see,  the  error  does  not  behave  monotonically.  (The  computations  above  were 
performed  by  PROBE,  which  offers  this  extraction  technique  feature).  An 
essential  prerequisite  of  the  extraction  here  is  the  knowledge  of  the 
coefficients  otj  and  |frj(0)  in  (2.8)  (and  the  adjoint  of  ^j(e)).  As  shown  in 
Section  2,  these  are  available.  For  the  extraction  of  other  data  of  interest, 
we  refer,  for  e.g. ,  to  [67].  Although  they  are  not  computationally  trivial, 
extraction  techniques  can  potentially  save  a  large  amount  of  computer  time 
(when  included  as  a  standard  feature  in  a  program),  especially  in  3  dimen¬ 
sions,  where  an  error  of  order  IX  in  stain  energy  is  easy  to  obtain  but  an 
error  of  order  IX  in  the  energy  norm  is  very  difficult  to  achieve. 
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8.  A-posterlorl  estimates. 


An  essential  aspect  of  the  finite  element  method  is  quality  control 
of  the  computed  results  of  interest.  The  p-version  (on  properly  designed 
meshes)  gives  an  effective  tool  for  this,  because  it  computes  extensions 
without  changing  the  mesh  (and  saves  precious  user  time).  The  computed 
sequence  of  the  value  of  interest  can  be  analyzed  by  various  extrapolation 
approaches  (or  simply  by  assessing  the  changes  with  p  visually).  If  the 
data  is  monotonic,  as  in  the  case  of  the  energy,  the  extrapolation  technique 
is  very  effective. 

Table  8.1  shows  the  approximate  relative  energy  norm  error  estimates  by 
the  program  PROBE  for  the  model  problem  mentioned  in  Section  3,  with  n  *  2 
layers. 


Table  8.1  The  estimated  and  true  energy  norm  errors. 


H 

N 

Estimated 

Error 

True 

Error 

P 

N 

Estimated 

Error 

True 

Error 

B 

41 

25.41 

25.41 

5 

497 

1.41 

1.47 

H 

119 

8.45 

8.45 

6 

695 

1. 15 

1.32 

B 

209 

3.91 

3.93 

729 

0.89 

0.98 

B 

335 

2.09 

2.  13 

n 

1199 

0.74 

0.85 

The  estimates  are  computed  by  using  extrapolation  based  on  the  formula 


(8.  1) 


cn'B-eex-efe 


where  (respectively  Epg)  is  the  exact  (respectively  computed)  finite 

element  energy.  (Note  that  there  are  3  unknowns  in  (8.  DiC.p.E^.  ).  We 
compute  out  of  three  successive  values.  The  final  value  of  Eg^  accepted 
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is  then  used  for  all  p.  When  the  error  curve  is  concave  (as  we  would  like  to 
achieve  by  the  proper  mesh  and  an  exponential  rate  of  convergence),  then  the 
estimated  errors  are  upper  bounds.  In  any  case,  as  Table  8.1  shows,  the  error 
estimate  is  of  high  quality.  We  cam  and  should  use  other  quality  controls, 
e.g.  various  equilibrium  checks,  etc.  (see  e.g.  (681).  PROBE  has  vairious  such 
features.  Essentially,  one  ha^  to  compute  the  values  of  interest  more 
accurately  than  needed  for  engineering  purposes  because  of  quality  control 
reasons. 

If  the  values  au-e  not  monotonic  than  it  is  easiest  to  present  the  entire 
sequence  to  the  user.  As  an  example,  we  show  in  Fig. 8.1  the  3  dimensional 
amalysis  of  a  splice  and  depict  the  maximal  principle  stress  in  the  region. 

The  standard  h- vers ion  computation  results  are  also  given.  The  data  are 
taken  from  [69]. 

It  is  observed  that  one  may  decide  by  inspection  that  the  p-version  has 
converged  satisfactorily.  A  similar  deduction  is  not  possible  with  the 
h-verslon. 


N(GlobalDOF) 

Fig.  8.1  The  accuracy  of  the  maximal  principal  stresses  in  a  splice 
computation 
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Similar  principles  can  be  used  in  the  error  assessment  of  the  modeling 
mentioned  in  Section  6  and  extraction  techniques  mentioned  in  Section  7. 

9.  Robustness  and  problems  Involving  "locking". 

A  robust  method  is  one  which  performs  uniformly  well  for  a  broad  class 
of  input  data.  Consider  for  example  the  elastic,. cy  problem  as  defined  in 
Section  2.  When  v— >1/2,  the  form  degenerates  (although  the  solution  stays 
smooth,  see  Section  2)  and  we  have  div  u— >0.  It  is  well  known  that  the 
h-version  of  the  finite  element  method  for  low  degree  p  performs  very 
poorly,  due  to  the  phenomenon  called  "locking". 

There  are  various  types  of  locking.  The  above  mentioned  is  called 

Poisson  locking  and  will  be  briefly  discussed  here.  Essentially,  when 

v  — ■»  1/2,  problems  arise  in  a)  the  convergence  in  the  energy  norm  and  b) 

computation  of  the  pressure  (o-  ♦  <r  ). 

x  y 

We  have  to  distinguish  between  the  case  of  straight  and  curved 
elements.  It  has  been  shorn  in  [70]  that  the  rate  of  convergence  in  the 
energy  norm  of  the  p-version  with  straight  triangles  is  not  influenced  by 
v  — > 1/2.  In  [71],  it  has  been  shown  that  the  h-verslon  with  straight 
triangles  does  not  show  locking  when  pk4.  The  general  case  is  analyzed  in 
[72],  [73].  Theorems  9.1  and  9.2  are  specialized  versions  of  the  more  general 
results  obtained  in  these  references. 

Theorem  9. 1  (3.3)  holds  for  straight -sided  triangles  and  parallelograms 

uniformly  in  v  provided  that  p 2 pQ  where  Pq  ■  4  for  triangles,  3  for 
Qp  elements  and  5  for  elements. 

The  above  theorem  shows  that  with  straight  triangles  and  parallelograms, 
no  locking  occurs  for  pkp^  both  with  the  h  and  the  p-verslon.  For 
curvilinear  elements,  we  have  for  the  p-version: 
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RELATIVE  ERROR  IN  ENERGY  NORM  Hell 


Theorem  9.2.  Let  u^  €  H^Q).  k>  1.  Then  the  following  estimate  holds  for 
the  p-version  uniformly  in  v, 

'Sfl-SFE^E  S  C<P"S>k  l|S0lHk(O) 

where  s>0  depends  upon  the  mappings  of  the  curvilinear  elements  onto  the 
standard  elements,  provided  these  mappings  are  rational. 

A  related  theorem  for  the  case  when  the  mappings  are  analytic  may  be 
found  In  [721,  [731. 

As  a  simple  illustration,  we  show  the  relative  error  in  the  energy  norm 
for  various  straight  and  curvilinear  choices  of  a  single  element  for  v  *  0.3 
and  v  -  0.5  -  10~1C\  where  a  type  element  is  used.  For  a  detailed 

analysis  and  numerical  examples,  we  refer  to  [721,  [731. 

3? 


NUMBER  OF  DEGREES  OF  FREEDOM  = 

Fig. 9.1  Error  behavior  for  curved  elements  for  v  *  0.3  and  v  *  0.5-10  *  . 

The  second  problem  is  the  pressure  recovery.  In  [741,  [75],  it  is  shown 
that  the  stress  components  <r^  -  <r^,  are  accurately  computable  directly 

from  the  solution  (by  Hooke’s  Law)  but  the  directly  computed  pressure  is 
unusable.  Nevertheless,  the  pressure  can  be  accurately  recovered  by  a  post¬ 
processing  technique  based  on  the  observation  that  it  is  a  harmonic  function 
and  that  Upg  is  accurate  in  the  energy  norm. 
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A  similar  behavior  occurs  for  other  cases  of  material  which  can  be 
associated  with  nearly  degenerate  bilinear  forms  (see  [14])  In  two  and 
three  dimensions.  The  h-p  version  Is  also  very  robust  In  relation  to  the 
shear  locking  that  occurs  in  plate  theory. 

10.  Implementations!  aspects. 

In  contrast  to  the  h-version,  the  p-version  needs  much  more  computational 
work  to  construct  the  local  stiffness  matrices  and  load  vectors.  Moreover,  it 
leads  to  matrices  which  are  less  sparse.  On  the  other  hand,  the  local 
stiffness  matrix  computation  is  completely  parallel  (and  for  uniform  p  is 
well  balanced)  which  can  obviously  be  implemented  by  parallel  computers. 

For  complex  geometries,  curvilinear  elements  with  relatively  large 
distortion  cannot  be  avoided.  This  problem  is  overcome  in  the  p-version  by 
using  quadrature  rules  with  the  number  of  quadrature  points  depending  on  the 
distortion. 

The  system  of  equations  for  the  FEM  solution  is  less  sparse  for  high  p 
than  for  low  p.  Hence  the  solution  is  more  expensive  for  high  p.  Never¬ 
theless,  the  ratio  of  the  computational  work  to  the  accuracy  obtained  is  more 
favorable  for  the  p-version  (this  is  also  true  for  engineering  accuracy).  For 
a  detailed  analysis  we  refer  to  [76],  [77], 

Iterative  method  techniques  like  the  conjugate  gradient  method  can  be 
very  favorably  Influenced  by  the  correct  selection  of  the  shape  functions. 

For  various  aspects  of  the  influence  of  the  shape  functions  on  the  iterative 
process,  we  refer  to  [78].  In  [79]  we  have  shown  that  in  2  dimensions,  the 
preconditioned  conjugate  gradient  method  (preconditioning  by  p  *  1  computa¬ 
tion)  requires  0(log  p)  steps  asymptotically  when  the  shape  functions  are 
properly  chosen. 
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In  (80],  detailed  experimentation  and  analysis  of  the  factors  influencing 
effectiveness  and  parallelization  on  Alllant  computers  are  presented.  It  is 
shown  that  the  speed-up  is  at  least  90%. 

In  general,  mesh  generation,  especially  in  3  dimensions,  is  a  difficult 
task.  Presently,  a  mesh  generator  geared  to  the  needs  of  the  p-version  is  not 
available  and  a  PATRAN  interface  is  usually  used. 

For  experimentation  with  mesh  refinement  of  the  h-p  version  on  tensor 
product  meshes  in  two  dimensions,  we  refer  to  [81]. 

11.  Engineering  experience  and  practice. 

In  the  previous  section,  we  discussed  various  theoretical  aspects  of  the 
h-p  version.  A  large  amount  of  engineering  and  industrial  experience  with 
the  method  has  been  gained  in  connection  with  the  use  of  commercial  programs 
FIESTA,  PROBE  and  research  program  STRIPE.  For  some  articles,  we  refer  to 
[82]  and  references  therein.  See  also  [83]. 
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under  the  general  administration  of  the  Director,  Institute  for  Physical 
Science  and  Technology.  It  has  the  following  goals; 
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